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Abstract 

The purpose of this work is to find the time dependent distributions of directions and posi- 
tions of a particle that undergoes multiple elastic scattering. The angular cross section is given 
and the scatterers are randomly placed. 

The distribution of directions is found. As for the second distribution we find an exact 
expression for isotropic cross section in 2 dimensions. For the same cross section we find its 
Fourier-Laplace transform in 3 dimensions. For the general case we devise a method to compute 
the Fourier-Laplace transform with arbitrary precision. 

This work is based not on the Boltzmann transport equation but on an integral equation 
formulation of the problem. The results are general in the sense that any initial condition is a 
linear combination of the cases considered in this article. 

Key words: Boltzmann transport equation, multiple scattering, integral equations, 
Lorentz gases. 
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1 Applications and overview 

As a quantitative analytical technique multiple scattering is used to study clouds [1], 
pollution, chemichals in solution [2], radiation in the human body [3,4,5], etc... A great 
effort has been made to study the scattering of laser beams by clouds. To this end 
some authors [1] have developped numerical stochastic models, based on the same general 
philosophy as ours, which have been shown to be equivalent to the transport equation 
which is more commonly used. Simulations are a major tool in the study of multiple 
scattering [1,6], but real situations are computationally so demanding that numerical 
methods based on analitical understanding of the problem can be necessary. 

In Physics multiple scattering is obviously related to Lorentz gases [7] and to diffusion 
of neutrons in nuclear reactors [8,9]. It has many applications in solid state [6]. In most 
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situations the particle being scattered is a photon or an electron. Here we shall work with 
a scalar particle, but the methods can be generalized for the case of internal degrees of 
freedom as long as interference effects remain negligible. 

In the 2nd section we find the directions taken by the particles of a beam which en- 
counters randomly distributed diffusing centers. Conversely the differential cross section 
of the particle-scatterer interaction can be solved in terms of the distribution of direc- 
tions. Asymptotics are also discussed. The results of the 2nd section are experimentally 
applicable as long as the detectors are placed at a distance to the medium to be studied 
which is much larger than the diameter of the medium. In the 3rd section the motion 
of the center of mass is found. In the 4th section we study the case of isotropic cross 
section. The probability density for the particles is found explicitely in 2 dimensions and 
in 3 dimensions its Fourier-Laplace transform is found. In the 5th section we show how to 
calculate the Fourier-Laplace transform of the probability density with arbitrary precision 
in the general anisotropic case for both 2 and 3 dimensions. 

The initial condition for all the results in this article is a beam of particles. This is the 
general case in the sense that any initial condition can be decomposed in terms of beams. 
The mathematical setting of this work is that of integral equations. For some instances 
this work provides an alternative to the Boltzmann equation for the study of transport 
processes. 

2 Diffusion of directions 

Consider the repeated elastic scattering of a point-like particle by randomly distributed 
recoilless diffusing centers. We assume that the scattering itself takes place in a region 
of negligible volume (or area) in the sense that the trajectory of the particle can be well 
approximated by a sequence of straight line segments. Let / be the angular cross section, 
where / is a function of the angle. If / is not equal to a constant, i. e., if the scattering 
is not isotropic, then a statistical description of the multiple scattering becomes difficult 
because the process is not Markovian in the physical space. However, the stochastic 
process for the directions is Markovian. In fact if, in d dimensions, we represent the 
direction along which the particle is moving by a point on S^ 1 , then the probability 
of transition from one point of S^ 1 to another will depend only on the angle between 
them, so that convolution techniques can be applied. We will show this in 2 and then in 

3 dimensions. 

In 2 dimensions let F(t, ip) be the probability density for the particle to be moving 
along the direction ip at time t. ~ will denote Fourier transform and Tp will be the Fourier 
frequency associated to the angle p. If A is the probability of collision per unit time then 
in an infinitesimal time dt the probability density F will change to: 

F(t + dt, p) = (1 - A dt)F(t, <p) + \dtj dip' F(t, p')f(p - </?'), (1) 
The Fourier transform of the above equation is: 

F(t+dt,Jp) = (l-Xdt)F(t,Jp) + Xdt F(t,ip)f(Jp) = F(t,Tp) + X dt F(t,p)(f(jp) — 1). (2) 
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This yields the differential equation: 



|lnF(^) = A[/>)-l], (3) 



whose solution is: 



F(t,ip) = F(0,<p) expA[/» - l]t. (4) 
An analogous reasoning shows that when A or / are time dependent: 

F(t,lp) = F(0,lp) exp fdt \{t) [f(t,lp) - 1]. (5) 

Jo 

These formulae are valid for any cross section /. In particular left-right symmetry 
(f(ip) = /(27r — ip)) is not required. 

When all particles move along the direction tp = at t = 0, F(0,Tp) = 1 and the 
above formulae simplify. This kind of initial condition, which we shall call stream initial 
condition, is important because an arbitrary initial condition can be constructed in terms 
of it. In 3 dimensions the stream initial condition will be a beam of particles moving 
towards the positive direction of the z axis, due to the definition of polar coordinates in 
3 dimensions. 

In 3 dimensions the equation corresponding to (1) is: 

F{t + dt, n) = (1 - A dt) F{t, VL) + \dt J dQ! F(t, Q') f(Q - Q'), (6) 

where Q — Q' is the angle between the directions Q and Q'. Here we are going to assume 
that the scattering probability / depends only on the polar angle 9, or, in other words, 
that if the particle undergoing scattering has internal (pseudo)vectorial degrees of freedom 
(such as the polarization of a photon), these are not involved in the scattering or are 
averaged over. We also assume, for simplicity's sake, that the initial condition F(0,Q) is 
independent of ip. Then F(t, Q) is also going to be ip independent and both F(t, Q) and 
f(Q) can be expanded in terms of Legendre polynomials Pi as follows (see [10] or [11]): 



00 21 + 1 

F(t,n) = Y,\l^ r F l (t)Y l \9) (7) 



1=0 

and 



/(n)=av /,li0(9) ' (8) 



1=0 

where 



Y?(e) = J^Pi(c<*e) (9) 



are the spherical harmonics and 



F,(f) =JdQ F(t, Q)Y l °(Q), (10) 
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and a similar formula for //. 

If we now substitute these expansions in the convolution product in (6) and use the 
sum formula for Legendre polynomials, 

and 

JdnYr*(n)Yp\n) = 8 w 6 mm ,, (12) 

we obtain: 

J dsn F(t,n') m - n') = E y ^ *i(f) /1 n°(n). (is) 

Proceeding as in the 2-dimensional case we obtain the time evolution of the coefficients 

my- 

F t (t) = F,(0) exp A(/, - / = 0, 1, 2, (14) 

as well as an equation analogous to (5). For the z-stream initial condition, Fi(0) = 1 
for I = 0,1,2,.... Different derivations of eq. (14) have been given by Goudsmit and 
Saunderson ([12], [13]) and by Lewis [14] (these derivations are also reviewed in [6]). 

Note that the fact that f(9) is arbitrary allows us to treat cross sections like the one of 
the Rayleigh scattering, which a Fokker-Planck equation on the sphere (e. g. the Debye 
rotational diffusion equation [15]) can not approximate. 

It is intuitive to think that asymptotically the directions will be uniformly distributed. 
From 

F(t,<p) 4E^.^) e X[f ~ m - 1]t (15) 

and 

/»= r&pf{<p)e-**, (16) 
Jo 

we see that unless f(ip) is a set of Dirac deltas placed at the vertices of a regular polygon, 
our final distribution in 2 dimensions is indeed going to be uniform. 

In 3 dimensions normalization implies fo = 1. We now show that the other coefficients, 
fi, 1 = 1,2, are smaller than 1: 

fi = JdQ y^^/(Q)^°(Q) = 2tt £ d(cosfl) f{9) Pi(cos6). (17) 

The Legendre polynomials take absolute values always smaller than 1 except for cos 9 = 
±1, where the even ones take the value 1 and the odd ones take the value ±1. Thus if 
f{9) 7^ a 6(9) + a 5(9 — n), where a + a = 1, then 
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/, < 2tt d(cos#) f(6) = Jdtt f(Q) = 1, 1 = 1,2, (18) 

and the time evolution (14) implies that in 3 dimensions the final distribution F(oo, Q) 
will be the uniform distribution. 

To measure quantitatively how the distribution approaches its final state, we choose 
the L 2 distance between the current distribution and the final one. In 2 dimensions this 
is: 

/. n^-d- (19) 

The Fourier transform is an isometry of the L 2 norm, so that the above distance is 
equal to: 



— V e 2X ^-^. (20) 

The time dependence of the second cumulant of the probability density for finding a 
particle which is suffering multiple scattering tends to be linear as t — > oo. By studying 
the approach to this behaviour Boguha et alii [16] found the above formula with the sum 
restricted to Tp = ±1. Higher terms in (20) are related to higher moments. 

In 3 dimensions the L 2 distance between the current distribution, F(t,Q), and the 
final one, can be expressed in a form analogous to (20) using expansion (7), the 
orthonormality relation (12) and the evolution equation (14). The final expression for the 
L 2 distance is: 

g^±l Fi(0 ) 2 e 2 ^- 1 )* (21) 

1=1 ^71" 



3 Motion of the center of mass 

If the cross section / is not ballistic then at every collision some of the initial momentum of 
the particle is transfered to the scatterer. As a result the velocity of the center of mass goes 
exponentially to zero. It is interesting to find how far the center of mass will go before (at 
t = oo) it stops. To study this note that in 2 dimensions for a particle (which we assume of 
unit mass) impinging along the direction f upon a scatterer, its momentum changes from 
(v cos (p, v sin if) to, on the average, (y / 27r dip' f (if' — (p)cosip',v J 27r dip' f \ip' — ip) sin if'). 
Thus Newton's law for the center of mass reads: 

f=Xv 
dt 

J df F(t, f) ( J df'f(f' — f)cosf' — cosf,J df'f(f' — f)smf' — smf\, (22) 
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or, using (15), 



r-2-K / (-271 (-271 \ 

y d(pe lflip lj dip'f((p' — <p) cos <p' — cos <p, y dp'f(p' — ip) simp' — simp \ = 



7T 7T 

/(^)vr[^_i + <^,i] - 7r[^_i + <5?,i], /(p)-[<^,-i - - -[^,-i - <%> i] 



(23) 



Since the initial condition F(0,t) and the scattering cross section /(</?) are real valued 
functions, F(0, -1) = F(0, 1)* and /(-l) = /(!)*. Therefore 



f=M^o,i)|^>-*(i,i) 



(24) 



For the initial condition p(0) = (i> Jq W dp F(0, <p) cosp, t> Jq W dp F(0, <p) sin the so- 
lution to this differential equation is: 



p(t) = Re[vF(0,l) e A[/(1) - 1] '(l,i) 
In particular for the stream initial condition and symmetric cross section, 

p(t) = v eW'W- 1 * I 
Finally, for x(0) = a last integration yields: 

x(t)=Re\vF(0,l) Ml,i) 

U L A[/(l) - 1] 1 ^ 

For the stream case with symmetric cross section this implies that the center of mass 
penetrates a distance 



(25) 



(26) 



(27) 



krloo 



(28) 



A[l-/(1)] 
into the scattering medium. 

In 3 dimensions for a particle of unit mass impinging along the direction Q = (9, ip) 
upon a scatterer, its momentum changes from (t>sin 9cosp>, v sin 9 sin p, v cos 9) to, on 
the average, (v £ dff sin 9' dp'f(ff -9,ip' - p) sin 9' cos p', v /* dff sin 9' dip'f(9' - 
9, p' - p) sin 9' sin p' , v £ dff sin ff ^ dip'f(ff -9, iff - p) cos ff). Thus Newton's law for 
the center of mass reads: 
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-1 = \ v r dOsinO r d<p F(t,9,ip) 
dt Jo Jo 

J dO' sin 9' J d(f'f(9' — 9, ip' — if) sin 9' cos <p' — sin 9 cos ip, 

PIT e-l'K 

\ d6 ; sin & \ dip / f(9 / — 9, tp f — <p) sin 9' sin cp f — sin 9 sin (p, 
Jo Jo 

j dff sin 9' j cV/(0'-0,¥/-p)cos0 / -cos0j. (29) 

As in section 2 we restrict ourselves to F and / functions of 9 only. Then only the 
z component of the above vector is non zero. Using the expression (14) for the time 
evolution of the Legendre polynomials coefficients of F(t, 9) we obtain: 



(0,0,/ dO'smff dtp' f (9' -9) cos 9' -cos 9) (30) 

JO JO 

Using the summation formula for Legendre polynomials (11) and the ortho normality 
of the spherical harmonics (12) we obtain 

p(t) = vF 1 (0)^~ 1 H. (31) 

For the z-stream initial condition, F(0, Q) = ^"^"^ , Fi(0) = 1, and p(t) = ve^'^k. 
From (31), 

2{t)=vF 1 (0) { ) k. (32) 

and 

KcoW-^jE (33) 
For the ^-stream case we obtain a formula similar to the 2-dimensional one: 

|i?(oo)l = A[l - A] 1 (34) 
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4 Isotropic scattering 



When the scattering is isotropic there is a direct way of obtaining the ^-stream probability 
density (p x ) in terms of the probability density for isotropic initial conditions (p): 

p x (t,x) = e~ xt 5{x - vt)5(y) + A /* dt' e~ xt ' pit - t', \x-vt'i\). (35) 

jo 

The rationale for the above equation is that at each collision the stream gives rise to 
an isotropic probability cloud. 

For the 2-dimensional case p is known [17]: 



p(t,x) = e 



-At 



5{r — vt) 



+ 



A 



2?rr 2nv J{vt) 2 - 



exp {j-\J (vt) 2 — r 2 ^H(vt — r) 



(36) 



where H is the Heaviside function. The integral in (35) can be done analitically (see the 
Appendix). The result for p x is: 



Px(t,x) 



-At 



8(x - vt)S(y) + 2M *_ x) exp {-J{vt) 2 -r 2 )H{vt 



r) 



(37) 



p x (t, x) is, in 2 dimensions, the probability density for finding the particle in x at time 
t provided that at time t = it left the origin moving to the right. As long as the particle 
doesn't collide its p x (t,x) will be 5(vt — x) 8(y), where v is the speed of the particle. 
Since the particle does suffer, on the average, A collisions per unit time, the density of 
probability contributed by the possibility of no scattering is e~ xt 5(vt — x) S(y)- Once the 
particle is scattered at a time V and to a direction ip' - which will happen with a probability 
density A/(c//)e -A *' - its contribution to p x (t,x) will be p x (t — t f , R- V '(x — vt'i)), where 
R-^i is the rotation that takes the tp' direction back to the positive x-axis direction and 
% is the unit vector that points to the right. All this paragraph is the following integral 
equation: 



rt r27r 

p x {t,x) = e~ xt 5(vt-x) 5{y) + \ / dt' e~ xt / dip' fUp') p x (t -t' , R-^ix - vt'i)) (38) 

Jo Jo 



The Fourier transform that we use for continuous variable is: 



/oo 
dx function(z) e~ l27VUX , (39) 
-oo 

and the same transformation except for the sign of the exponent for its inverse. The 
Fourier transform of the integral equation is: 

p x (t, V) = e -^ e ~W)^ + A /* dt' r dip' fUp') e- xt 'e- i2 *W* vt 'p~ x (t - t' , R^P), (40) 

JO JO 
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where {V) x stands for the x component of V. We now take the Laplace transform (denoted 
by a hat) and obtain: 

/■27T 

(u + A + i2-k(v) x v) p x (u, V) = 1 + A / dip' f(ip') p x (u, R-^V), (41) 

Jo 

where to is the Laplace conjugate variable to time. 

Consider the case of isotropic scattering (f(p') = Then the integral in the preced- 
ing equation is the Fourier-Laplace transform of the solution to the problem of isotropic 
scattering with isotropic initial conditions. For this case the following integral equation 
can be written for its probability density p [8]: 

p(t, x) = rj(t, x)+\J* dt' J dx rj(t', x) p(t -t',x- x) (42) 

where 

V(t,x) = ^5(\x\-vt). (43) 

Using the Fourier and Laplace convolution theorems one can solve for the Fourier- 
Laplace transform of p [8] : 



P = (44) 
1 — Xrj 



in terms of the Fourier-Laplace transform of rj: 



f}(u, u) = 1 (45) 

^(A + uoY + (2ttw) 2 



Substitution into (24) yields: 

1 



Px{oj,v) 



(uj + A + i27r(P) x v) (1 — \r](uj, V)) 



(to + A + i2n(v) x v) + uf + (2ttH 2 - ^ ' 

which is the Fourier-Laplace transform of (37). 

In 3 dimensions the analogous procedure starts from the equation: 

p z (t,x) =e~ xt 5{x)5{y)5{vt-z)+\ J* d£ e~ xt ' J dQ' f(9') p z (t-t' , R^(x-vt'k)). (47) 

We are assuming here that the scattering is independent of the azimuthal angle (p. R-ci' 
is the rotation which maps the Q' direction into the z-direction and has its axis in the 
xy-plane. The scattering cross section is normalized as follows: 
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J dQ f(6) = 2nJ*dO sin0 f{9) = 1. (48) 
The Fourier-Laplace transform of the integral equation is now: 

(u + A + i27v(u) z v) /5 z (u, V) = 1 + A J dQ'f(9') ^(co, R- a ,v), (49) 

Consider the case of isotropic scattering (/(#) = ^). Then the integral in the preced- 
ing equation is the Fourier-Laplace transform of the solution to the problem of isotropic 
scattering with isotropic initial conditions. For the this case the integral equation (42) 
holds and its solution is again given by (44). The difference is that r) is now: 



and fj(u, v) is: 



e -xt 

V(t,r) = — 5(r-vt), (50) 



v) = - arc tg— — . 51 

1-KVV A + ui 



Finally, the Fourier-Laplace transform of p z (u, v) in 3 dimensions for the case of isotropic 
scattering is: 

Pz[U ' V) = (u + A + i2n(v) x v) 27rw-Aarc tggf ' { } 

5 Anisotropic scattering 

The form of the integral in (41) is similar to a convolution. This suggests the following 
procedure to factorize the unknown function out of the integral. First, write V as R^V: 

/■27T 

(u + A + i2ir(R (p v) x v) p x (uj, R^V) = 1 + A / dip' f(ip') p x (u, R^-^V). (53) 

J 

Define now the following integral transform for a function g of x (x is a 2-dimensional 
vector): 

g{x,Jp)= dip g(R^x) e~^, JpeZ. (54) 

J 

From Fourier series we know that this transform can be inverted: 

g(R v % = ±- *£g(x,lp)e**>. (55) 

271 TpeZ 

It has the following property: 

g(R a x,ip) = g(x,ip) (56) 
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Thus the transform that we have defined is, up to a phase factor, the Fourier transform 
of the restriction of the function to a circumference about the origin. 

Application of Jq W dip to both sides of the integral equation (53) yields: 



— — r Z7r ~ ._ 

27r5o,^ + Xf(Tp)p x (u, v,Jp) = (cu + A) p x (u, V,Tp) + ilrtv j dip p x (u, R v v) {R v v) x e~ %w '. 

J 

(57) 

Let (pa denote the angle between z7 and the x axis. Then the integral in the above equation 
can be developped as follows: 



/ dip p x (uj,R v V) (R^V) X e 

J 

/ dip pJ^.RJu) \V\ 
Jo 



uptp 



i(<P+<Pa) _|_ p-i(<P+V>a) 



e**P x (u,a,Tp- 1) + e- i ^p x (cu,u,ip+ 1) 



(5f 



Substitution of this result back into (57) yields a set of recursive equations for the trans- 
form p x of the probability density: 



27T<J (W + [A(/(y) - 1) - u] p x (u, V, ip) 



VKV\V\ 



e^p x (u,i7,ip--l) + e-^p x (u;,v,ip+l) 
If the scattering has left-right symmetry (f(ip) = f(2n — ip)), then, 

/» = /(-?) 

For /5 X this symmetry implies: 
From this and (56): 

p x (w,u,±<p) = e ±l ^ p x {u, \v\T,<p). 
Now the set of equations (59) can be simplified to: 

2ir-up x (u, | v\,0) = i2-nv\V\ p x (u, \v\i, 1) 

and 



ipeZ. 



(59) 



(60) 



(61) 



(62) 



(63) 



Wiv) - 1) - /5a.(w, ¥>) = i7Tu|i?| p x (w, |z?|?, ^ - 1) + p x (u, \P\l, ip+1) 



for (p=l,2,..., 



(64) 
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or (bearing (61) in mind and remembering that /(0) = 1) simply to 



M 1 ~ /(¥>)) + w] Pxiu, <p) + i™\V\ \p x (uj, (p-l) + P x (oj, \P\i, up + 1) = 2irS ^, 



for = 0,1,2,.... 

This more compact form brings out the matricial structure of the equations: 



(65) 



i7rv|z/| [A(l — /(2)) + oj\ mv\V\ 

mv\V\ [A(l — /(l)) + oj\ mv\V\ 

2mv\v\ [A(l - f(0)) + u] J 



p x (w, |z?|?,2) 
p x (w, 1) 

V p x (w, i^r, o) j 



(66) 



/ • \ 




V 2tt J 

From (55), (61) and (62) the Fourier-Laplace transform of the probability density is: 



p x (u,u) = -\ + 2^ p x (u, \u\l,ip) COS((fi (fff) 

7T \ Z — in 



(67) 



Thus by approximating equation (66) by a finite dimensional truncation we can find the 
Fourier-Laplace transform of the probability density with arbitrary precision. There are 
recursive formulae that allow to invert finite tridiagonal matrices at a low computational 
cost [18]. 

We assume that at t — all particles leave from the same point and have the same 
speed. If this is not the case, summation over speeds and positions is necessary. Let 
I if) = -^zZipeZ Hp) e ^ be the initial angular distribution of velocities. Let p (u,v,Tp) 

be the function corresponding to p x (u>, v,Tp) when the initial direction of the particle is ip. 
R-p is the adjoint operator of R^, therefore: 

Then, using standard properties of Fourier sums and property (56) of the ~~ transform, 

p2tt 



— — 

J o dip I((p) p v (u, V,Jp) = I(Tp) p x (u, u,Tp). 



(69) 



In 3 dimensions we can proceed analogously. The starting equation is now (49). In it 
we write V as R^V. 



{u + A + i2n{R n ») z v) p z {u, R n P) = 1 + A J dtt'f(6') p~ (u;, R n -n>v), 



(70) 
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3-dimensional rotations can not be added commutatively, but Rn-w is well defined. It is 
the rotation that maps the Q' direction into the Q direction and has its axis perpendicular 
to the plane defined by these directions and the origin. 

Define now the following integral transform for a function g of x (x is a 3-dimensional 
vector): 

g(x,l) = J dVl g(R n x) Y^ty, UN. (71) 

If g does not depend on the azimuthal angle p the following 2 properties hold: 

a) 



g(R Q xJ)=g(xJ)J^- [ Y l (Q) (72) 



and b) the above transform is one-to-one and its inverse is: 



g(x) = ^°*(0, 0)=E\p7^ m 0, (73) 

UN leN " ^ 

where (0,0) means (6,(p) = (0,0). Similarly to (10) we define: 



Application of / gJQ ^°(^) to both sides of the integral equation (70) and use of the 
convolution theorem (13) yields: 

(u + X) p z (ou,v,l) + i2nv J dttp z {uj,Rnv) {R^) z Y?(Q) = v 7 !^ + A/ Z p 2 (cu, u, I). (75) 
If Q = (0, (/?) and the polar coordinates of V are (9p, if a), then geometry yields: 

{Rq,v) z = \v\ {cosQcosQp — sva.9sva.9pcos{ip — (pp)). (76) 

Due to the assumed cilindrical symmetry of the scattering cross section we can take 
(pt; = 0. Then, writing the trigonometric functions as linear combinations of spherical 
harmonics, we obtain: 

{R,V) Z = ^ \P\ [l?(n)cosfc- (Y^Q) -y^n))^]. (77) 

When we substitute this expression back into (75) we are going to have products of 
spherical harmonics. By means of Clebsch-Gordan coefficients [19] the products Y^ zl (Q) YJ°(fi) 
can be written as linear combinations of ^±i(^)- But it easy to check that for a cilin- 
drically symmetric function g(x), J dtt g(Rnx) Y l m (Q) = if m ^ 0. Therefore only 
the product ^°(^) remains. Using the recurrence relations between Legendre 

polynomials [10] this product can be expanded as follows: 
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v /(2/ + l)(2/ + 3) 



v /(2/ + l)(2/-l) 



(78) 



Substituting this result back into (75) we obtain the following set of recursive equations 
for p z : 



FW1 „, , - , _ n i27rdz7|(Z + 1) -£ ^ i27rdz?U n ~ , -> , lX 

[A(l-//)+wJ p 2 (u;, z/, Z)H — = cos 0? p 2 (a;, z/, Z+1)H — . = cos 6*^ p z (u, u, l-l) 



(2/ + l)(2/ + 3) 



^/(2Z + 1)(2Z-1) 



In matricial form: 



for 1 = 0,1,2,.... 



(79) 



^3 cos Qv [X (i-f 2 ) +u ] 








i2wv\u\2 
VIE 



COS 0j7 



^pcos^ 
[A(l-/i)+w] 

COS 0* 







COS fj7 



i2irv\i?\ 

[a(i-/ )+w] y 



V p z (^,o) y 



/ • \ 
o 





(80) 



where, for symmetry's sake, we have not simplified [A(l — fo) + oj] to oj. 

Truncation of the equations to a finite / yields a set of equations that can be solved. As 
in the 2-dimensional case the matrix is tridiagonal and it can be inverted using recursive 
formulae [18]. By means of the inversion formula (73) we can thus find the Fourier-Laplace 
transform of the probability density with arbitrary precision. 

As in the 2-dimensional case we assume that at t — all particles leave from the 
same point and have the same speed. Let = Y J i, m h,mXi n (Si) be the initial angular 

distribution of velocities. Let Pq(oj, P, I) be the function corresponding to p~ z (w, P, I) when 
the initial direction of the particle is fl. R-n is the adjoint operator of Rq, therefore: 



Psi(v,v,l) = p z (u,R-nP,l). 



(81) 



Then, using standard properties of spherical harmonics and property (72) of the trans- 
form, 



J dn i (n) p n (u, p,i) = J p z (u, v, i) ii, 



(82) 
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This result does not depend on the value of the lj )m 's for m 7^ 0. As noted earlier 

/ dVL I(Q)p n (ui), V, I) can be inverted by means of (73) to yield the Fourier-Laplace trans- 
form of p(t, x) only when I is cilidrically symmetric. Thus to apply formula (82) one has 
to decompose a general initial condition as a sum of cilindrically symmetric (with respect 
to different axis) components and then apply (82) to each component separately. 
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Appendix 

Substitution of (36) into (35) yields 3 terms, which are: 

e- xt 5(x-vt)5(y), (83) 

Xe -*l , df W-«®-«l-<» (84) 
Jo 2ir\x - vt'i\ 

and 



r t exp(±J(v(t-t')) 2 -\x-vt'i\ 2 ) 
A e / dt' J-H(y{t - f) - \x-vt'i\). (85) 

Jo 2nv^/(v(t-t')) 2 - \x-vt'd 2 



Since the integration variable is not the argument of the Dirac delta function in the 
second term, the following formula will be needed: 

/&«(«(*))= E w55T (86) 



u(x)=0 



It is convenient to write the integral in the second term (84) as follows: 

J_r> SUZ-ti-vt + v)) 
2irv Jo 2n\x — yi\ 

The derivative that appears in the r.h.s. of (86) is then: 

-?-(\x-yi\-vt + y) = l+ l~ X (88) 
oy \x — yi\ 

Using of the last 3 formulae the second term becomes 

Ae~ A * 1 



2-kv y — x + \x — yi\ 
The argument of the Dirac delta function in the second term vanishes for 



(89) 
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Substituting this in (88) yields 



Denning 



Xe' xt 1 
2nv vt — x 



A = vH 2 - r 2 



and 



B = 2v(x-vt), 
the third term (85) can be written as 



2irv 
Xe~ xt 



vt-r oo 

r^y dt , 

Jo t^n 



3-1 



X 2 e~ xt rU^k ,., ^X j (A + B_t'Y- 



2nv(x — vt) 
Xe- Xt 



XV A + Bt> ' 
cxp 1 



2(vt-x) 



2nv(vt — x) 
Collecting results we obtain (37). 



fXVvH 2 -r 2 \ 
exp - 1 



16 



References 

[I] A. V. Starkov, M. Noormohammadian, U. G. Oppel, "A stochastic model 

and a variance-reduction Monte-Carlo method for the calculation of light 
transport", Applied Physics B, Lasers and Optics, B 60, (1995) 335-340. 

[2] Milton Kerker, "The scattering of light" (Academis Press, 1969). 

[3] R. F. Bonner and R. Nossal, Appl. Opt. 20 (1981) 2097. 

[4] A. P. Sheperd and P. A. Oberg, Laser-Doppler Blood Flowmetry (Kluwer, 
New York, 1989) 

[5] B. Chance et alii, Pro. Natl. Acad. Sci. USA 85 (1988) 4971. 

[6] J. M. Fernandez- Varea et alii, Nucl. Instr. and Meth. in Phys. Res. B73 
(1993), 447-473. 

[7] R. Garcia-Pelayo, Physica A, 216, (1995) 299-315. 

[8] G. I. Bells and S. Glasstone, "Nuclear Reactor Theory" (Van Nostrand 
Reinhold, 1970) 

[9] J. H. Ferziger et al, "The theory of neutron slowing down in nuclear reac- 
tors" (Pergamon Press, New York, 1966). 

[10] Harry Hochstadt, "The Functions of Mathematical Physics", Dover (1986). 

[II] Chun Wa Wong, "Introduction to Mathematical Physics", Oxford (1991). 
[12] S. Goudsmit and J. L. Saunderson, Phys. Rev. 57 (1940) 24. 

[13] S. Goudsmit and J. L. Saunderson, Phys. Rev. 58 (1940) 36. 
[14] H. W. Lewis, Phys. Rev. 78 (1950) 526. 

[15] C. S. Johnson Jr., Don A. Gabriel, "Laser light scattering", Dover, (1995) 
p. 53. 

[16] M. Boguha et alii, Physica A 230 (1996) 149-155. 

[17] J. Masoliver, J. M. Porra and G. H. Weiss, Physica A 193 (1993) 469-482. 

[18] Y. Huang and W. F. McColl, J. Phys. A (1997) 7919-7933. 

[19] C. Cohen-Tannoudji, B. Diu and F. Laloe, "Quantum Mechanics", Wiley- 
Interscience (1977). 



17 



